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1. Introduction 


This work represents part of a project to develop an atmospheric general circulation model based on 
the semi-Lagrangian advection of potential vorticity (PV), with divergence as the companion 
prognostic variable. The resulting model derived by Bates and Li is hence termed the PV-D model. The 
motivation for the use of potential vorticity as a prognostic variable is that it is the most basic 
dynamical quantity (see [5]) and is conserved along a trajectory in frictionless adiabatic flow. It is 
reasonable to assume that, particularly in situations of highly nonlinear flow, use of the evolution 
equation for PV as one of the governing equations will give a more accurate simulation of the PV 
evolution than a model using velocity components or vorticity and divergence as the basic prognostic 
variables. An important consideration in the development and time discretization of the PV-D model 
was the idea that a multigrid (MG) solver could be developed for efficient solution of nonlinear implicit 
equations obtained at each time step. The focus of this part of the project was the development and 
testing of the multigrid solver. 

In this report, the continuous and discrete equations will be presented, ending with the system to be 
solved by multigrid at each time step. The development of the MG solver will be outlined, the 
algorithm described in detail, and results on model problems will be presented. 

2. The Model Equations 

Here, the equations will be presented. Much of the detail and motivation for choosing the model is 
outside the scope of this report and will not be described in detail. In the model, a generalized vertical 
coordinate £ is used, rather than an isentropic vertical coordinate. This hybrid coordinate basically is 
related to the ratio of pressure to surface pressure in the lower atmosphere (and thus is constant on the 
surface) while it is isentropic (follows layers of constant temperature) in the upper atmosphere. The 
reason for this is to eliminate problems at the lower boundary layers where stratification may become 
nearly neutral. Closely related to the above, potential vorticity is replaced by a quantity referred to as 
pseudo potential vorticity (PPV) and it is this quantity that is advected. In the upper layers, where the 
vertical coordinate is isentropic, this quantity is basically potential vorticity and is conserved. This 
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quantity is not exactly conserved in the lower layers. However, it can be argued that the isentropic PV 
does not have much dynamical significance there. Physics is also included in the equations as moisture. 

2.1 The Continuous Equations. 

A longitude-latitude coordinate system is used in combination with the vertical coordinate, giving a 
system (X, <j>, 4). A € [0,2k] is longitude, <j> e [-k/2,k/2] is latitude (with <j> = -k/2 as the south pole 
and <j> = n/2 as the north pole), and £ is the vertical coordinate, with bounds depending on the 
definition chosea There are eight primary unknown functions defined on this domain: 

p Pressure. 

if/ Streamfunction. 

X Velocity Potential. 

M Montgomery Streamfunction. 

9' Potential Temperature perturbation. (Potential temperature 0is written as 9+0 , here 
#(4) is a given mean potential temperature and 6\X, <f>, Q is the perturbation.) 

4 Vertical velocity. 

o) Pseudo potential vorticity. 

q Specific humidity. 

Physical constants that appear in the equations are defined as follows: 

a Radius of the earth (6.371229-10 6 m). 

Q Rotation rate of the earth (7.292- 1 O' 5 1 /s). 

C K 0.28704/1.005 (non-dimensional). 

C p 1.005-10 3 J/kg-K. 

poo 10 5 Pascal (N/m 2 ). 

R Universal gas constant (287 J (°K)'' Kg' 1 . 

k RJCp. 
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Other variable quantities are: 


n Defined as C p {plp<mf and used in the equations for convenience. 

/ Coriolis force (2 Q sin$. 

P \!a df/dtf> = 20 cos <pla. 

The quantity (3 is related to the semi-Lagrangian discretization and will be discussed fiirther in the next 
section. As noted below, the “continuous” form of the divergence equation given here is already 
discretized in time, and includes /?. 


In the following, d/dt denotes the Lagrangian, or material, derivative (along a particle path). V is the 
horizontal gradient operator 


1 d- \_d-_ 
a cos <j> dX' a d<f> 


where a is the radius of the sphere. V 2 is the horizontal Laplacian: 


( 2 . 1 ) 


V 2 =- 


1 d 2 


a 2 cos 2 (j) d?i a 2 dtp 2 

Although the more familiar quantities of velocity and divergence are not explicitly included in the 
equations here, they are important and it is useful to see how they relate to the unknowns used. A 
Helmholtz decomposition of horizontal velocity is used to obtain the streamfonction ^and velocity 
potential x •> s° that u and v satisfy 


(2.2) 


u 1 d% \ dy 

acos<j>dX a d<f> ’ 

\ dx \ dy/ 

v = — — + — . 

a d<j> a cos <f> dX 

The horizontal divergence D satisfies 

D = V-(u,v) t = V 2 x- 


(2.3) 


(2.4) 


Now, the set of continuous equations considered here can be presented. 
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The PPV equation 


da d$ 
— = a — — 
dt d$ 



1 d$ d 

( \ dx 1 dxy' 

1 d'$d_ 

( 1 dx 1 dy/\ 

a cos <j> dX di J 

yd d<f> a cos <p d, l j 

a d<t>d£ 

^acosfdX a d<f> ) 


a cos^ 


<?TI <3(1 + 0.608 q)0 <?TI^(l + 0.608g)fl 
dX d<j> d<f> dX 


0p!H) 

Here, S\ is a term due to turbulent friction. 


(2.5) 


The Divergence Equation 

The divergence equation is given here discretized in time. This is because the equation is derived by 
first performing a semi-Lagrangian time discretization of the horizontal momentum equations, then 
taking the divergence. An uncentering parameter is used and x\ depends on this parameter and the 
timestep At. Superscript n+\ denotes values at the new time step. 

(v^r + r,(v 2 A/)" +1 - h u n+l (v 2 e) n+ ' - r,(vnr‘ .(v«r 

-(vn)" +1 + r 1 0.608(6h) n+1 -(v 2 !!)" 


\/i+l 


+ 1 


l &x l dy 


\/7+l 


^ r '' a cos <f> dX a d<j> 
- forcing 2 = r 2 + S 2 At. 


-A i( v V)” +1 


( 2 . 6 ) 


The Hydrostatic Equation 


dM 

T-,dO n , nn 

R ' 

J 

p ri 


-n — + 0.608 


q6 

— i 

d$ 


\ Poo ) 


v Poo J 


dp_ 

d$ 


= 0 . 


(2.7) 


The Continuity Equation 


d_ 

dt 


dp 

{#) 


r dpY 


v 2 *+ 


d? 


= o. 


( 2 . 8 ) 
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The Thermodynamic Equation: 


dO a Qi 

Tr Q '=7v*^ 


(2.9) 


The Definition of f 


Z(p,p s ,ff) = t. 


( 2 . 10 ) 


The Definition of ox 


(£)-- 


/ + VV 

dpjdl; ' 


( 2 . 11 ) 


The Moisture Equation: 


d? __C& 

dt L ’ 


( 2 . 12 ) 


where Qi is the apparent moisture sink and L is the latent heat of condensation for water vapor. 


Boundary Conditions. The boundary conditions needed to close the system are as follows: 

= = (2.13) 

M(t o ) = mj) + (0(Zo) + 0(Zo)) n(&), (2.14) 

and 

= ( 2 - 15 ) 

Here, pr is a given function, generally taken as constant in tests here. 


3. The SLSI Discretization 


The form of the equations used (particularly the divergence equation) is tied to the use of a semi- 
Lagrangian vector discretization. The idea of semi-Lagrangian discretization is to express the material 
derivative of a function in terms of the difference of the function's value at a point at the current time 
step (the "arrival" point) and at the corresponding "departure" point at the previous time step, obtained 
by backtracking along the trajectory. In practice, the departure points cannot be determined exactly, 
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since these depend on the new, undetermined velocities, and exact integration along trajectories may 
not be possible. Instead, the departure points are approximated by assuming that the average velocity 
along a trajectory is its value at the midpoint, and obtain these midpoint values by spatial interpolation 
of time-extrapolated velocities from previous timesteps. To make this more precise, suppose that a 
time step A / has been chosen, and the solution is known at time nAt. A two time-level semi-implicit 
scheme is used, giving equations relating the solutions at times nAt and («+l)A/. The time at which a 
quantity is defined is denoted by superscripts. Now, let x be a point at time level n+1 , and let x> denote 
the corresponding departure point at time level n. Using the trajectoiy midpoint approximation for 
velocity V gives x*=x-At V^ l2 ((x+x-)/2), where velocities at time step n+l/2 are obtained by forward 
extrapolation in time, with V t+l7 =1.5V’-0.5V This implicitly defines x*. Tthe superscript n + 1 denotes 
function values at arrival points at time level n+ 1, while ( .)♦ denotes the function values at the 
corresponding departure point at time level n. When the spatial discretization is introduced, the level 
H+l values will be defined at the grid points, while the departure point values will be obtained by 
interpolation to the departure point at time level n. 


A straight forward semi-Lagrangian discretization of the divergence equation leads to Rossby wave 
instability (see [2]). For this reason, as in previous work, instead of using the value of the Coriolis 
parameter / at the departure point at time step n, we use a Taylor series expansion of f 


= /(0-A0vA /+/' 


(3-1) 


where P(<f>) = (\! a) df [dip, v is the ^-component of velocity approximated at the trajectory midpoint. 


and /' represents the higher-order terms in the expansion. (Numerical experiments show that it is 
sufficient to keep only the second-order term in /’.) v is approximated in the first-order term by the 
average (v" +l + v" )/2 , and in / byv.” . 


7 



The Grids: Here, the computational and the physical grids are introduced. As noted previously, the 
computational domain is given by [0,2rc]x[-7t/2, tt/ 2] x[<^,£max], where <$> and depend on the 
definition of the. vertical coordinate chosen. Several examples for the definition of 4 are given in 
Section 3. Vertical discretization is based on K levels, which are prescribed so that 4b < 4i < — <44 ~ 

^„ax . In the following, these are referred to as the “solid” levels (for historical reasons), and can be 
irregularly spaced. Intermediate levels (“dashed” levels) are defined as 4* = (4*-i + 4* ) 1 2 for k = 

1 For placement of one of the equations, “virtual” levels are also defined by 4* = (4* +^ k+i )/2 
for k= \,2,...JC-\, and 4 0 = 4o* Note that for uniform spacing, the virtual levels and solid levels 
coincide. 

Now, for each fixed 4-level used, the horizontal discretization is defined. A rectangular discretization of 
the (T,4) plane is chosen as follows. Let Nx> 0 and N,>> 0 be fixed integers, which define the mesh 
spacing on this grid by hx= 2n/Nx and h#= n/N ^ The computational grid points are then defined to be 
(A;, 4, 4), where X, = i-h\ for i= 0, l,...,Nx and 4 = -n/2+j- h# for j= 0, 1, . . ., A^. (At times, fractional 
indices are used, and the meaning should be clear. For example, 4 + 1/2 — (4 +4+i )/2.) I n the present 
work, Nt = N is taken as some power of 2, and Nx — 2 N, resulting in a uniform horizontal 
computational mesh with hx~ h^=h = n/N. In the physical domain, each 4-level is a sphere of radius a 
(the increment due to height is neglected). Since X and (j> represent longitude and latitude, respectively, 
(A 0 , 4 , 4 ) and (/f?Ar, 4 ,£) correspond to the same physical point for each j, and, for all i, ( T/,4? , 4) 
corresponds to the south pole while (X, ,4v, 4) corresponds to the north pole. Thus, on the sphere, the 
corresponding grid has uniform latitudinal (north/south) mesh spacing of size a ■ h, while the 
longitudinal (east/west) spacing depends on the latitude and is given by a h cos 4> so that it equals the 
north/south spacing at the equator and approaches zero towards the poles. (Generally, it is not 
necessary to distinguish between the computational and physical grids, but keep in mind that the 
equations are defined on the sphere.) 

Positioning of the Unknowns and Equations. The unknowns and equations are all defined along the 
same vertical lines (that is, there is no horizontal staggering of the grids used for each unknown). 
However, the vertical positioning differs from unknown to unknown and equation to equation. The 
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equation positioning will be specified later as the discrete equations are presented. The solid levels 
(§’s) are used for defining p, 6 (and thus, <9 and 0), £, and q. The dashed (intermediate) levels are 
used for defining M, x, y, and ox One exception is that M 0 is used to denote the value of M at the 
surface (a solid level) although it is defined by a boundary condition there, rather than an equatioa For 
those functions defined at solid levels, subscripts ij, k correspond to the value at point (A.,, $, &), while 
for those defined at dashed levels, these subscripts indicate the function value at ( A i , , 4 ) • 

The Discrete Operators. First, the discrete form of the vertical derivatives used, along with some 
necessary notation, is introduced. An approximation to the vertical derivative 0- /dQ is used primarily in 
equations centered at the dashed levels. Letting A be defined at the solid levels, this approximation is 
defined by: 

fa) AzAl. ( 3 . 2 ) 

bk bk-\ 

A vertical derivative is also used at the virtual levels in (3.14) below. This time letting A be defined at 
the solid levels, this approximation is defined by: 

Ism) (3.3) 

Finally, because Mo is defined at the surface (a solid level) and M\ is defined at the first dashed level \ , 
the special placement of equation (3. 14a) makes the following definition necessary: 

{~8m) (3.4) 

At dash ed levels, vertical averaging of some quantities at solid levels is also needed, and is denoted and 
defined by 

2,={U + 4-,)- (3- 5 > 

The quantities denoted by (.) k are defined at the virtual levels ~% k by quadratic interpolation from solid 
levels as follows: 

O* = »i.., (•)*-, + ^,oU +»i,,(.)* +1 (3- 6 ) 
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where 


(6 -&]("&-&«) 
*'° (6 -&-.)(& -5k + .) 


(3.7) 


Note that with equal vertical spacing, Wk,o — 1 and Wt ,. i — Wk,\ — 0. 


The horizontal derivatives are defined as in earlier work on shallow water equations (see [8]). These 
are fairly straightforward second-order finite differences away from the poles, and the integral 
definitions of the operators (where possible) at the poles. In the following, keep in mind that domain is 
periodic in the ^.-direction, so that a subscript of ij,k is the same as i±2NJ,k. Now, the first-order 
differences away from poles (that is, where 0 <j < N) are defined by: 



A+ij^k A-\j,k 

2a cos tj>h 


(3.8) 


and 


M,,, - 


A. 


j+\,k iJjJzl iL 

lah 


(3.9) 


Note that, although A and <j> are used in the notation, these are derivatives on the sphere, and not in the 
(A, 0) plane. For the horizontal Laplacian V 2 , a conservative discretization is used in order to avoid 
problems over long time integrations. For the sake of simplicity (although at the risk of confusion) the 
same notation is used for both the continuous and discrete operator. Away from the poles, we define 


_ Ai+iJJk ~ ^ Ai,j,k Ai-\j,k 
' a 2 h 2 cos 2 4>j 


COS ^y + j/2 ( y+l,A ■ Ai,j,k) ~ cos < Pj-l/2^.Ai,j,k ' Aij-\,k ) 
a 2 h 2 cos <j)j 


(3.10) 
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At the poles, an approximation for the integral definition of V 2 is used. Since the first index is 
meaningless there, it is omitted in the pole value itself. The operator at the north pole (j = N) is defined 
as follows, with a similar definition for the south pole: 





(3.11) 


In the equations, some quantities appear in the equations that involve derivatives with respect to X and 
0, which are not defined at the poles, although these can be well defined and exist in the limit. For such 
quantities, an average over points adjacent to the pole is used. Now, the discrete equations can be 
presented. For simplicity, the subscripts ij have been omitted, and it should be understood that these 
equations hold (with the appropriate discrete operators) for all ij and the indicated values of k. The 
forcing terms, labelled forcing 1 , ... ,forcing6 , are computed based on quantities at the previous 


timesteps. 


Pseudo Potential Vorticity Equation: This equation, the discretization of (2.6), is defined at the dashed 
levels, and holds for k= 1 ,2,-X Note that equation (2. 1 1) has been used to eliminate co from the (2.6), 
so that pseudo potential vorticity no longer explicitly appears in the equations at all. This simplifies the 
system somewhat. Here, q\ and qi arise from semi-implicit treatment of the pseudo-potential vorticity 
term, and are given in terms of the solution at previous timesteps. It should be noted that <72 = 0 at the 
poles. 

[i-zf(^](/+vV*)-Tf[^(n)^((i+o. 608 9 x^+^))] t 

+ *f[^(n)^((l + 0.608? )(0+ 3 12 

+ $5 x S-VERT; - ifS^-VERTj 

+[(?,)* +(q 2 )k(^<Pk +^z k )](S 4 p) =(, forcings k . 
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Divergence Equation: The discretization of (2.7) is also defined at the dashed levels, and holds for 

w 2 Xk + zfv 2 M* - <[nv 2 e\ - 

S l O.6O%[s x nS x (q(0+ 0)) + Sfl-S+(q(0+ (3-13) 

if 0.608[<?( 0+ 0)V 2 n] k +/3 t i(s xXk ~ S+<p k )-fS l V 2 <p k = ( forcing 2) k . 


Hydrostatic Equation (surface) : This equation, the discretization of (2.8) near the surface, is centered 
at (0.75^0+0.25^,): 




— fg + — T, + ^1 ( )\ + a 2 {d/:0)2 


+0.608/d 


3i - n 

J q(0+0)- i 
H P) o 


+ Mq(0+0) 


n 


[afd^p\ + a 2 (S 4 p) 2 ] = (forcingl) 0 . 


(3.14a) 


Hydrostatic Equation: This equation, the discretization of (2.8) away from the surface, is centered at 
the virtual levels and is defined for k— 1 1 : 


(s i M) k -|n*[r*+r* +1 +(<%*)* +(<%*) 4+I ] 


+— 0.608 
2 


( - 


n 


K^q(&+ 0)— 


[( S 4 p) k +(^p) k+ 1 ] = {forcing!) k . 

k 


(3.14b) 


Continuity Equation: This equation, the discretization of (2.9), is centered at the dashed levels and is 
defined for for 

1 + f(v 2 Zk) + = (forcing*)*- ( 3 - 15 ) 

Thermodynamic Equation: This equation, the discretization of (2.10), is centered at the solid levels 
and is defined for for k=0,\,2,...,K. 

{0) k + ^f k ^ k =(forcing5) k . (3.16) 
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Definition of vertical coordinate: This equation, the discretization of (2.1 1), is centered at the solid 
levels and is defined for for kr=\,2,...JC-\: 

&Pk’Po’0+ ff k) = £k- ( 3A? ) 

The Moisture Equation: This equation, the discretization of (2. 1 3), is defined for for k= 0, 1 ,2 
Since this explicitly defines q at the new timestep, q can be taken as a known function in the remainder 
of the equations and this equation can be ignored in the development of the solver. In all tests, <7 = 0. 

q k = (forcing^) k . (3.18) 


In the above equations, a number of quantities have been used, which are defined here. 

if A + fiVfe 


a i =■ 


=- 




(3.19) 


6 - 6 

The VERT and VERT terms in (3.12) are approximations to S$u and <%v, respectively, and are defined 
as: 


and 


VERT x u = a^8 f (s xX -S, ¥ )\ +b^{s xX ~ S^\, 

(3.22) 

VERT f = a k \s,(s a ~ + **[M 8#- 

(3.23) 

VERT? = a K [s { (s^-S^)] Ki + b K [s,(s xZ - S.y)]^ 

(3.24) 

VERT? = a,[^(^+ S x yr)\ + 6 ¥ )\, 

(3.25) 

VERT ; = «,[*<(**+ 

(3.26) 


(3.27) 


where the coefficients a and b are interpolation coefficients to the proper levels as follows: 


6 - 6 . 


a, = =H — = 


6,-6 


b — 6 _. 

6-6 
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(3.28) 


a - & §k_ - 

u k ~ - ’ 


b k = & I*" 1 , k = 2,3,..., AT — 1; 

_ £*-i 


% 


_^/r ^y-2 . ^ _ _£at ‘tat-i 

£a:-i ^a:-2 ^k-2 ~ £k- i 


Boundary Conditions: The boundary conditions for the discrete equations are: 


ii 

n 

© 

(3.29) 

m q -< t>uj)+(e 0 + # 0 )n 0 . 

(3.30) 

Pk ~ Pt(^"> $)• 

(3.31) 


In summary, the system under consideration consists of the equations (3. 1 2) - (3. 1 7), with boundary 
conditions (3.29) - (3.31), and the set of unknowns being ip, O', y/, %, M). In tests, two definitions 

for £ in (3.17) are used. The first, called the cr- coordinate, can be used for testing with with no 
orography, and is given by: 

Z=o{p,p 0 ,p r ) = ]- . (3.32) 

Po~Pt 

The second, a hybrid coordinate, is defined in terms of cr above, and is given by: 

1 - eT aa 

£=3nin + S(^)-[0--3mn] where g(£T) = - — — . (3.33) 

1 6 

Here, 6Ln is some value less than the maximum value of 6, and in tests is set to 265. a is a value that 
controls the transition between the surface fitting ^-coordinate and the upper isentropic layers, and in 
tests is taken variously as 5, 10, and 20. The values for rare basically half the time step At, possibly 
modified by a mild uncentering parameter. In all tests, the fs are taken as 1 800 s. This completes the 
definition of the discrete equations. 

4. The Multigrid Algorithm 

In this section, the multigrid (MG) algorithm for the system given in the previous section is developed 
and described. It is assumed that the reader is fairly familiar with multigrid methods, and is referred to 
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[3], [4], [6], [7] and [9] for more detailed descriptions of some of the processes. It should be noted that 
the general philosophical approach taken here is to err initially on the side of robustness, with lull 
efficiency as a secondary consideration. 

In Section 4. 1 , some general terminology, notation, and algorithms are introduced, avoiding as much as 
possible the details of the multigrid components. Most of these details are given, along with the reasons 
for the choices, in Section 4.2. The choice of relaxation is quite involved. Some general decisions are 
given in Section 4.2, but the main task of analyzing the system and deciding the final form of relaxation 
to be used are contained in Section 4.3. Results contained in Section 5 show that for the given 
discretization, a reduced grid version of the algorithm may be needed. Details of this version, in which 
the fewer grid points are used along latitude lines near the poles, are given in Section 4.4. Finally, some 
comments on improving the algorithm, particularly with respect to efficiency, are given in Section 4.5. 

4.1 Terminology and Notation 

To describe the MG algorithm, it is helpfiil to first introduce some general terminology and notation. 
The terms “grid” and “level” will be used interchangeably. Superscripts will generally be used for the 
level number, with the finest level indicated by 1, and the coarsest by M. The level / grid is denoted by 
Cr, and consists of the points with 0<i < 0< j < N‘ t , and 0< k < K 1 . When the level is 

understood, the point ( X\ , ^ , g ' k ) is also referred to as point (i,j,k). Because of the diverse vertical 
positioning of the various equations and unknowns, it is more accurate to say that the level / problem is 
defined relative to this grid rather than on this grid, but this should cause no confusion. The coarser 
grids will be defined later. Now, let u denote a solution approximation on level /. Here, u l consists of 
the approximations (or corrections) to all the individual grid functions p, 9, %,%, y/, and M, defined at 
the appropriate points. Similarly,/ denotes a level / right-hand side, corresponding to the different 
equations (3.12) - (3.17), and again defined at the appropriate points. Interpolation (prolongation) from 
level /+1 to level l is denoted by P , while restriction from level / to level /+1 is denoted by r! + x . P 1 is 
used for interpolating corrections from coarser to finer levels, while R* +l is used for transferring 
residuals to coarser levels. In the FAS algorithms (see Section 4.2), a different restriction (generally 
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injection) may be used for transferring the values of the unknown funetoions to coarser and is denoted 
by R l+ ' . In the FMG process, a higher-order interpolation (here, cubic) is used for interpolating the 
approximation to a new finer level. This is denoted here by P 1 . Now, let F* denote the discrete 
approximation to the full nonlinear operator on level /, so that the fine grid problem can be written as: 

F'(u') = f'. (4.1) 


In the MG algorithm, it will be necessary at times to look at grid functions and operators defined on 
only part of the grid in order to perform line or plane relaxatioa Subscripts will be used to denote 
partial grid functions and operators. A line (ij) in Cr is the vertical line of points (ij, k) with 0 < k < K 1 . 
This is also denoted by G/ ■ . A y'-plane is the vertical plane of points (ij,k) with 0 < i < N[ and 

0<k<K‘, and is denoted by Gj . Let u/ y denote the part of u defined the vertical line (ij), and u' } be 
the part of u defined on the y'-plane. For completeness, let d N be the part of the function lying on the 
vertical line at the north pole, and u' s be the part at the south pole. Separating the operators is 
somewhat more complicated. For this purpose, the complement of a partial grid function is introduced. 
These are denoted by a in front of the subscript, and indicates the remainder of the grid function. 

For example, consider the example of the y-plane (The following notation also applies to an (ij) line). 
The part of u 1 defined everywhere except on plane j is denoted by , so that u can be written as 

(Uj,u[j). Now, let F l . denote the part of r defined on the plane j. This is actually defined on the full 
function u, but the primaiy interest is on the action of F‘ on u' r Thus, Fj(u'j) should be interpreted 
more accurately as F^u'^u 1 ^ ) with frozen. For the plane relaxation, some linear approximation of 
Fj around the current approximation z/ y is needed. This is denoted by L y , so that for a perturbation 
Su’j of u‘j, iljdu’j » Fj(Uj +Su i j \u[ j ). Here, of course, L y is generally a function of u. For convenience, 
setting this linear approximation will be denoted, somewhat inaccurately, by Z, y » Fj (u( ; u( ] ) . In the 
current problem, the details of defining ll } will be given later. For plane relaxation, the problem to be 
solved for the correction « y to u( is then 

Wy = f* = fj ~ F j ( u j » u ~j) ■ (4.2) 


Similarly, (i,j) line relaxation requires the solution of a problem of the form: 


(4.3) 
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for a correction u‘ t J to u ‘ Uj , where as above, l! t J is a linear approximation to F'j ( u‘ t J , u[ t J ) around u‘ . 
Generally, a direct solver is used to solve (4.3). 


For solution of the linear y'-plane problem (4.2), a multigrid solver is used. The levels used here will be 
indicated by a second superscript, with the levels being denoted by /, /+1, ... , M. Note that the finest 
level number is the same as the global level number. This is somewhat less messy to write than starting 
again with 1 , and is actually more convenient in programming. In the work here, the same coarsest 
level number Mis used for the global and plane MG solvers. Let the local grids be denoted by 
Gj J ,Gj l+l . Because of the coarsening chosen (see Sections 4.2 and 4.3), interpolation and 

restriction (written as matrices) between levels m and m+ 1 depend only on m, and not on / or j. Thus, 
these can be denoted by P™ and /?" +1 , respectively. A sequence of linear operators iif for levels m = /, 

/+1, ... , M, are defined, with lif = l!j and the coarser grid operator defined as some standard 
approximation to Z/ y . When j and / are understood, they are dropped from the notation. 

4.2 Major Multigrid Choices 

Now, there are a number of decisions to be made in constructing the multigrid algorithm. Although the 
choices are interrelated, they can often proceed from the most general (overall type of algorithm) to the 
specific (form of relaxation). In this section, the choices are outlined and explained. As noted, we start 
with the most general choice of algorithm type. 

Full Approximation Scheme (FAS). To handle the nonlinearity, the lull approximation scheme (FAS) 
[3] is used, rather than manipulating the problem to produce linear equations using some outer 
linearization scheme like Newton's method. In FAS the nonlinearity is incorporated into the cycle and 
taken to all coarser levels. In this way, the full nonlinear problem can be solved in essentially the same 
time it would take to solve one of the linearized problems. This requires some reformulation of the 
multigrid algorithm, both in relaxation and the form of the coarse grid correction. Rather than 
formulate a coarse grid equation for the correction (that is, instead of using the residual equation), the 
coarser grid problems approximate the full fine grid equation, with an appropriate correction added to 
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the right-hand side. The type of relaxation needed is determined by examination of the linearized 
problem, and then modified for the nonlinear case. In the notation given in Section 4. 1 , ( 1 , 1 ) V-cycles 
using FAS can be written recursively as follows. The relaxation used will be defined more precisely 
later in this section. 

FAS Multigrid Cycle: FASMG ( ) 

If / = M Relax n c times on F l (u') = f l (call GR{M\u M \f M ): See Section 4.3). 

Otherwise: 

1. Relax on F'(u‘) = f 1 (call GR(l-,u';f')). 

2. Set u M = R ,+l u‘ and f M = R ,+ \f‘ - F‘(u‘)) + F M (u M ). 

3. Call FASMG (l + \;u M ;f M ). 

4. Set u 1 <-u l + P'(u M -R M u‘). 

5. Relax on F'{u‘) = f‘ (call GR(l;u' ;/')). 

Full Multigrid (FMG). For efficiency, it was decided to employ full multigrid [3], [4], [9], which is a 
process in which the first approximation on each level is obtained by interpolating the solution from the 
next coarser level (here, using cubic interpolation). In FMG, only a fixed reduction in the error (and 
therefore a small number of cycles) is needed on each succesively finer level to obtain a solution 
accurate to the level of discretization error, given that the coarser level solution is sufficiently accurate. 
Although the time discretization is separated from the solution of the nonlinear implicit equations, the 
time-dependent nature of the problem is important here. The FMG process is designed to solve the 
given implicit equations at each time step to the level of spatial truncation error. Since the 
discretization error also depends on At, the problem must be formulated to solve for the change in the 
solution over the previous timestep, rather than the full solution. Otherwise, error accumulation much 
larger than the discretization error can result. Solving for the incremental solution change is easily 
accomplished by using the solution from the previous timestep as an initial approximation, and by 
constructing the FMG process so that, on a given level, the initial approximation to this time increment 
is obtained by interpolation from a coarser level. This is done in the algorithm below. 
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Another consideration arising from the time dependent nature of the problem is the accuracy of the 
solution of the implicit equations that is required at each time step. It is easy to argue that, for an 
isolated equation or system, there is no need to solve below the level of discretization error since, 
beyond that point, the work invested does not result in a better approximation to the solution of the 
continuous problem. In time dependent problems, however, there may be accuracy or stability concerns 
because of error accumulation or possible feeding of nonlinear instabilities. Accuracy is addressed by 
the use of the FMG algorithm (designed to solve for the increment). This should ensure that overall 
error remains below the level of discretization error. Also, in this model, extrapolation of nonlinear 
terms, which is the (main) source of nonlinear instability, is avoided, so that no special care is needed to 
exclude or filter problem frequencies. Now, the FMG algorithm, based on vFASMG cycles per level 
and written for the increment, follows: 

Full Multigrid A Igorithm: FMG ( u 1 ; /' 1 ) 

1 . For / = 2, 3, ... , M, set u M = R M u‘ and f M = R l+l (f l -F l {u‘)) + F l+X 

2. Set l-M. 

3. Call FASMG {l\u‘ -J l ) v times. 

4. If / = 1, stop. 

Otherwise, set u‘ <r~u‘ + P‘(u M - R l+l u l ), ! <- 7-1, and go to 3. 

Global Grid Coarsening. The type of grid coarsening to use is a basic decision, here based primarily 
on experience and on knowledge of (and partly on ignorance of) the geometry of the mesh. A fairly 
obvious choice is to use standard coarsening in the A and (f> directions. Such coarsening is typical in 
spherical models. Coarsening in £ is more problematic, since little is known a priori about the strength 
of connections in the vertical direction relative to that in the horizontal directions, particularly since it is 
anticipated that vertical mesh spacing can be nonuniform, with higher resolution near the surface and 
top boundaries. Since horizontal coarsening alone reduces the number if grid points by a factor of four, 
relaxation work in a V-cycle is 4/3 times that on the finest grid. Coarsening in the vertical direction also 
would only reduce that factor to 8/7, while possibly requiring much more complicated (and expensive) 
relaxation schemes to get the smoothing required. Thus, it was decided not to coarsen at all in the 
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vertical direction. In the notation introduced earlier, this means that, letting N 1 = 2 J+l N, the level / grid 
is defined as follows: 

G , = {(A‘j l J ,4l)J = 0X...aN I J = 0X...,N l ,k = 0,\,...,K}, 
where A- = i-x/N l , tflj = -x/2 + j-n/N 1 , and £ = 4*. Note that (A' +l , <j> 1 ^ , 4* +1 ) = (A' 2i , <jx lj , ) for / = 1, 

2, ..., A/-1 . In the following, superscripts on £ values will be dropped, since the same vertical 
gridpoints are used on all levels. In tests, it was found that coarsening to N l = 4 was sufficient. Some 
effects of this grid on relaxation will be discussed in the next section. 

Interpolation and Restriction. Given the semicoarsened mesh, fairly standard interpolation and 
restriction operators can be defined. Since there is no vertical coarsening and horizontal coarsening is 
uniform, simple bilinear interpolation for each unknown can be used within horizontal layers. Similarly, 
horizontal foil weighting is used for restriction for the residuals of each equation. 

4.3 Relaxation. 

A final decision on relaxation details must be based on the discrete operators involved, since a basic 
MG principle is that relaxation must smooth the error in the coarsening directions. As will be verified 
later, the horizontal Laplacian (V 2 in the equations) is a principle part of the operator and is involved in 
smoothing. Some decisions can be made based on this knowledge, the geometry, and the coarsening 
used. When point-type coarsening is used, smoothing generally occurs in the directions) of strong 
connections (large coefficients). When the problem is coarsened in a direction of weak connections, it 
is necessary to simultaneously relax blocks of points that are strongly connected to each other (e.g., 
line relaxation). Smoothing then depends on the size of the out-of-block connections, which are now 
relatively large. This is important here because of the uniform coarsening in the latitude-longitude 
directions. Ignoring the vertical connections for the moment, in the neighborhood of the poles, a 
uniform latitude-longitude grid gives much smaller mesh spacing in the A (east-west) direction than the 
(f> (north-south) direction. Thus, for the V 2 operator, near the poles the connections in the A direction 
are much stronger than those in the <f> direction. For this reason, some form of coupled relaxation along 
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such latitude lines is necessary to obtain full horizontal smoothing of the error. The question of 
coupling of the equations in relaxation will be addressed later. 

As noted previously, little is known about the strength of vertical connections relative to the horizontal 
connections that may occur in general circumstances. To allow for strong vertical coupling (as well as 
to forestall possible problems due to vertical central differencing and averaging) a robust approach in 
which relaxation couples the unknowns on each vertical line is chosen. This coupling also influenced 
the choice to eliminate the pseudo potential vorticity ryfrom (2.5) using the definition of <y(2.1 1), 
which only quantities within a vertical line. 

Taken together, these choices indicate a multigrid approach in which (//) vertical line relaxation is used 
away from the poles (and at the poles themselves), while a 7 -plane relaxation is used near the poles. 
These will be defined later, after analysis of the system. 


Relaxation Analysis. In a complicated system such as this, it can be difficult to identify the principle 
part of the operator, which is the set of terms that actually determine the smoothing. Typically, this is 
done by writing the linearized system in operator form with coefficients frozen around the current 
solution approxiamtion, taking the determinant, and finding the dominant terms (see [4] and [10]). The 
operators contributing to these terms are used in relaxation (with other terms frozen, either for the 
entire sweep or just part of the sweep, such as relaxation on one line). For a problem such as this, the 
dominant terms can change according to the solution approximation itself (since the problem is 
nonlinear), the location within the domain, and the grid level, making a full analysis difficult. Several 
approaches were used to help in the identification of these terms, including systematic trial and error. 


A technique was developed to facilitate identification of these dominant terms. For this, the Jacobian is 
computed locally using the given solution approximation. To see how this works, the discrete system 
of n nonlinear equations can be written in a general form as follows (temporarily deviating from the 
notation introduced earlier by letting superscripts enclosed in parentheses indicate equation or 
unknown numbers): 


F </> (m ( 1 ) ,« ( 2 ) ,...,w < ' ,) ) = 0 (/ = l, 2 ,...,/i). 


(4.4) 
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where the u flJ (the /’ th unknown) is a grid function defined on the set of points Q / (identified by their 
indices (i,j,k)) and (the /’th equation operator) is a vector valued function defined at the set of 
points Q. 1 '. Denote the components of by and of£^ by The sets of points can be 

different for each unknown and equation, and, in general, a point identified by (ij,k) in one set need not 
necessarily coincide with the point with the same index in another set. However, all points and 
equations with the same ij lie on the same vertical line, and those with the same k are close to each 
other (within half a vertical layer). Now, we are interested in the linearized operator around (/ 'J'k). 
Since all discrete operators involved are local, a neighborhood A can be found such that each 

F?]. k . only involves unknowns defined in this neighborhood. Thus, for each l,m pair, a spatial array of 


values 


4 '& = 




(4.5) 


st- 
em be found (where the derivative is computed numerically, using calls to an existing residual 
calculation routine). This can be interpreted as the stencil of u <m) in F/, j, k . . For convenience, the 

subscripts are now dropped unless needed for . We are interested in how to smooth using the local 
operator, which we write symbolically as: 

V (U) . . . J (, - n) ' 

J = l . . . • (4.6) 

j(n, 1 ) 

As stated previously, the terms responsible for smoothing are those that dominate the determinant of J. 
Rather than actually multiply the stencils involved, a faster (and rougher) idea of the dominant terms in 
the determinant can be obtained by simply looking at the sizes of the operators involved. For this, we 
define 




max 

O',/',*) 




(4.7) 


the maximum stencil entry, and examine the resulting matrix 
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(4.8) 






The determinant of this matrix is taken, term by term, and those products with the largest values are 
noted. Here, only a few dominate, while most are zero. From this list, the important f ’ m> ' s are 
identified. However, each such stencil can have more than one contributing term in the equation Fy = 

0. Generally, these can be easily identified by matching the form of the stencil with the known form of 
the discrete operators involved (that is, it is easy to tell whether the stencil corresponds to a horontal 
Laplacian, or a vertical first derivative, for example). This gives a good idea of the main terms in the 
discrete operator responsible for smoothing in the neighborhood of ( i'J'k '). In addition, it can give 
information on whether or not some nonlinearities are important in relaxation, and whether coupling of 
equations is necessary. A note of caution is required here. At times, cancellation of what appear to be 
dominant operators can occur. This is not revealed in this analysis, so further examination of the 
equations is warranted if there are several large, nearly equal terms in the determinant. When such 
cancellation occurs, it may be important to couple relaxation of these equations/unknowns, and to also 
take into account the next largest determinant terms in relaxation. 

When the dominant terms involve vertical derivatives, this process needs further refinement. It is 
necessary to take the terms involved into account in relaxation, but this also means that other terms 
may be responsible for horizontal smoothing. Thus, the stencils are summed in the vertical direction 
and the process is repeated. This has the effect of eliminating those terms containing vertical derivatives 
(although not vertical averaging). The remaining dominant terms terms must then be identified and also 
included for relaxatioa This is somewhat more difficult, since the determinant of the corresponding 
system is generally zero, so large terms are identified by inspection. 

Such a process was quite helpful in the current work, although in practice, there was quite a bit of 
testing and verification to ensure that smoothing is adequate. Another useful technique was to compare 
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the residual of the relaxation system (corresponding to a line or plane) with the full residual there when 
the global solution is corrected. If the full operator residual stagnates while the relaxation problem is 
solved well, this indicates that some important terms are being neglected. 


For convenience here, equations (3.12)-(3.17) will be referred to as equations 1-6, and the unknowns 
(p. O', 4, y/, z M ) as unknowns 1 -6. Now, we can begin to analyze the smoothing behavior of a 
model system with no orography and a representative solution approximation. The mesh is defined by 
= 64, JV* = 32, and K = 1 8. A sigma vertical coordinate system (3.32) is used with uniform vertical 
mesh spacing (i.e., £ = k/K). The matrices (4.8) are computed at several places in the domain. For 
illustration purposes, the point (32,17,2) is first examined. This lies on the equator near the surface. 


Now | J | is given in the following matrix: 


2.34E-10 
2. 99E-12 
1 . 38E-01 
1.80E+01 
0. 00E+00 
1.01E-05 


4 . 34E-02 
0.00E+00 
0.00E+00 
3.21E+09 
8 . 32E+04 
0.00E+00 


2.79E-09 
9.10E-06 
8 . 75E+03 
0.00E+00 
1 . 00E+00 
0.00E+00 


3.29E-14 
1.02E-11 
0 . 00E+00 
1.82E-03 
O.OOE+OO 
0.00E+00 


1.02E-11 

3.29E-14 

0.00E+00 

0.00E+00 

O.OOE+OO 

0.00E+00 


0.00E+00 
1. 84E-08 
1 . 80E+01 
0.00E+00 
O.OOE+OO 
O.OOE+OO 


The determinant has six nonzero terms, given in descending order in the following table. 


Product Terms 

Size 

1,5 

2,4 

3,6 

4,2 

5,3 

6,1 

6.076E-17 

1,5 

2,3 

3,6 

4,4 

5,2 

6,1 

2.554E-18 

1,5 

2,6 

3,3 

4,4 

5,2 

6,1 

2.51 IE-18 

1,4 

2,5 

3,6 

4,2 

5,3 

6,1 

6.327E-22 

1,2 

2,5 

3,6 

4,4 

5,3 

6,1 

4.729E-22 

1,3 

2,5 

3,6 

4,4 

5,2 

6,1 

2.532E-24 


Each row of the table corresponds to one nonzero term of the determinant, where each term of the 
determinant is the product of one entry from row 1 , one from row 2, etc. The pairs of numbers in the 
table are the row (equation) and column (unknown) numbers of the factors in the determinant term, 
while the “size” column is the term itself (with sign omitted). Thus, for example, the first row 
corresponds to the term 




This is the largest term. Examination of the actual stencils in the Jacobian shows that these terms 
correspond to the operators shown below (with coefficients omitted): 

y (1 - 5) ~ v 2 

j< 2 - 4 ) ~ v 2 

J iX6) ~ s 4 
J (4 ' 2) ~ s e 

_ J 

Z 6 ' » ~ i 

This identification process is repeated at a number of locations throughout the domain. At this stage of 
the work, rather than to include the different terms only where they tend to dominate, the choice has 
been to attempt to ensure robust smoothing and include all terms that dominate anywhere in the 
domain. Once this preliminary examination was performed and candidate terms for inclusion were 
identified, tests were done to ensure that line and plane relaxation did, in fact, work as expected. This 
was done by comparing the full nonlinear residuals for the lines and planes before and after relaxation 
to see that they were sufficiently reduced. 

In testing plane relaxation, an FASMG solver was implemented, allowing for some nonlinearities. 
However, this was found to be unnecessary, and only a linear approximation to the equations was 
necessary. This linear systems used in line and plane relaxation can now be specified. The same terms 
are used in both types of relaxation, in all parts of the domain, and on all levels. The linearized 
equations are as follows, with the perturbations from the current solution designated by a single 
underline, and frozen coefficients designated by a double underline. All fully-frozen frozen terms are 
taken to the right-hand side, which is now denoted by / . 


1-rf (5f'<p k v V* - <(/ + v V + (9i)* + 


(#eP)k=f\ (4-9) 


V 2 z t + rfv 2 M k - ^[nv 2 ^ -Mv 2 <p k =f 2 


(4.10) 
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fe«) 0 -(fa + jn,][a,(^), + a ,(t s g.h] 


4 r ' +ai( ^ }i + a ^D^y Po h + \ n / Pl p^ 

(4.11a) 

40.608a: | ?(#+£)% j +J 9(0+£)%j [«i(^p), + a 2 (^p) 2 ] = /o 3 


(s ( u) k -|n k [(s 4 0 ) k +(s 4 ff) k+l ] 


-f [r* +T M HS ( ff) k +(s 4 ff) M ][[u/ P ]p ^ 

(4.11b) 

+i 0.608 A^(0 4 £)fV/^ [( 3{P)k + (^P)k+\ ] = fk 


4%p) 'J v2 *J + 4^p ) k (^) k +[i + <( v2 ^) + =/* 4 

(4.12) 

(i) k +Al=f k \ 

(4.13) 


(4.14) 


For line or plane relaxation, the V 2 operator is restricted to that line or plane. Now, given a line or 
plane, the linear operator L corresponding to the linearized F^ or Fj can be constructed. Thus, 
relaxation at one line (ij) on level / can now be defined: 

Nonlinear Line Solver: NLS(l,i,j,u',f') 

1. Define L » F i l j (ul J ;u l _ i J ). 

2. Solve Lu=f = f‘j - 

3. Set u' j <- u \ j + u . 

A direct linear solver is used for step 2. The efficiency of this step will be discussed later. 
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Plane Relaxation. For test purposes, a direct solver for the plane problems based on cyclic reduction 
was initially implemented. Not surprisingly, it proved to be very expensive, and therefore impractical. 
However, exact solution of the plane problems is not necessary, since the goal of relaxation is simply to 
reduce oscillatory off-plane error components. This can be done much more efficiently using a plane 
MG solver. Some problems were encountered here, and will be described in the next section. Consider 
the plane j on global level l. Using the notation introduced in Section 4. 1 , the grid is denoted by G \ , 
and consists of the points ( , </>] , & ) , with 0</<2A/' and 0< £ < A\ As inthe global MG solver, this 

grid is not coarsened in the vertical direction. Standard coarsening is used in the A-direction, giving a 
2D semicoarsened grid. As before, denote this fine grid by level / and denote the grids used in this 
plane solver by G‘j m , with m = l , ... , M Thus, G' m consists of the points (A?, #},&)» with 0</ <2N m 
and 0 < k < K. Note that, due to the coarsening used and the grid numbering chosen, level m uses the 
same values for JC[ and AT for the global grids and all plane grids, independent of /. Because of the 
semicoarsening, vertical line relaxation is used in the smoothing on each level, which is done in 
red/black order. Linear interpolation is used in the T-direction only. Similarly, A-direction full weighting 
is used for residual transfer to coarser levels. The coarsest grid problem (generally four vertical lines) is 
solved by several relaxation sweeps. The operators L m for m = /,..., M are defined here using the 
differences given in Section 3 in the linearized equations in (4.9)-(4.14). In the next section, an 
alternative is presented. Now, the plane MG (1,1) V-cycle has a form independent of / and /, and is 
defined as follows: 

The Plane Multigrid V-cycle: PMG(rn,u m ,f m , L M ) 

If m = M: Relax n c times on L m u m = f m (vertical line relaxation, red/black ordering). 

Otherwise: 

1. Relax on L m u m - f m (vertical line relaxation, red/black ordering). 

2. Set w m+l =0 and / m+l = R m+ '(f m - L m u m ) . 

3. Call PMG{rn + \,u m+l ,f m+ ',L m *\...,L M ). 

4. Set « m u m + PpU m+ ' . 

5. Relax on L m u m - f m (vertical line relaxation, red/black ordering). 
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Plane relaxation is defined in terms of the plane MG cycle above, and is given here. 

Plane Relaxation: PR(l,j,u' ,f'). 

1. Define iij * Fj(aJ;uij), and coarser grid operators . 

2. Set u‘= 0 and /' = // - F' («}; «i,). 

3. Call PMG(l,u‘ ,f‘ ,lif 

4. Set w' <- + u / . 

It should be noted that the linear systems used in line relaxation of the global sweeps and in the line 
relaxation used within the plane solver are basically identical. As noted, these line problems are solved 
directly, and the work there constitutes the major part of the overall solution algorithm. Efficiency 
considerations will be discussed further in Section 4.5. 

At the poles, relaxation uses vertical line solves very similar to the (ij) line solves away from the poles. 
The only real difference is in the discretizations used in the equations (4.9)-(4.14). For convenience, 
pole relaxation will be denoted using the same notation as the nonlinear line solver, with (ij) ~ (1>®) 
the south pole and (1,/V') for the north pole. Now, the full global relaxation sweep can be defined in 
terms of^ane, a parameter that specifies where line and plane relaxation are to be performed: 

Global Relaxation Sweep: GR(l;u l ;f') 

L Call NLS(l,\,0,u',f') and NLS(l,\,N l ,u' ,f') (pole relaxation). 

2. For j= 2, 4, ..., V- 2, 1, 3, ..., V- 1, 

IfKI^pi^call PR(l,j ,u‘ ,f‘) . (plane relaxation) 

3. For i+j odd, then i+j even (0 < / < 2 N 1 , 1 < j < N 1 - 1 ), 

If \4\ < ^ plane , call NLS(l,i,j,u',f‘). (line relaxation) 

This completes the specification of the basic MG solver. 
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4.4 The Reduced Grid. 


Because of problems with divergence using the current discretization discussed in Section 5, it is 
sometimes necessary to work with reduced grids. The basic idea is that the small mesh size in the X- 
direction near the poles is really an artifact of the uniform longitude-latitude grid, and does not 
contribute to overall accuracy. Thus, it is possible to seek a solution whose smoothness along latitude 
lines (from one grid point to the next) depends on the latitude line j. To make this more precise, 
consider the fine grid G l , with mesh parameter V . As before, the grid consists of the points 


with 0 < i < 2N\ 0 < j < N\ and 0 < k < K x . For this type of reduced grid, for each j, a number p is 
chosen, generally a power of 2, so that the X values actually used along line j are X,, i — 0, f>j , 2 p , ... , 

2 V. Let Nj = N'/Pj ■ Function values at intermediate points can be defined (if needed) by interpolation 


in the /l-directioa In the work here, linear interpolation was used. This affects the discretizations for 
the gradient and horizontal Laplacian operators given in (3.8) and (3. 10) as follows: 








2 a p } cos <f> h 


(4.15) 


and 


(v 2 ^l) = - ^!z£>±L 

' ‘J* a f^jh 2 cos 2 <j>j 

COS <p J+v2 { A,, j^,k ~ A,j,k) " C0S fij-I/liXij'k " Ai,j-\,k) 
a 2 h 2 cos <f> } 


(4.16) 


Note that this decreases the size of the connections in the ^-direction, and by choosing appropriate p's, 
it may be possible to eliminate the need for plane relaxatioa However, the approach taken here has 
been only to take the smallest p needed for convergence and retain the plane relaxatioa The actual 
choice of the p's will be discussed further in Section 5. Here, it will be assumed that p = 1 except at a 
few /-planes near the north and south poles. 


So far, only the finest global reduced grid has been defined. The method of coarsening such grids must 
be specified. For this, a superscript is added to indicate level number. That is, define p) ~ Pj. Now, the 
coarser grid p's can be defined recursively. Remember that plane j (with/ even) on level / corresponds 
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to the plane jtl on level /+1, the next coarser level. Let & indicate the reduced level / grid. Rather than 
“fully’’ coarsening each grid (so that the number of grid points along all even latitude lines, including 
those already reduced, is halved on coarser levels), we take G ,+l =& n G /+l . That is, if a plane/ on 
level / contains fewer than 2 V points in the A-direction, the number of points is not reduced further on 
going to the next coarser level. This can also be interpreted as p' + ' = max(l,/^ ; / 2) (or, alternatively, 

that Nj +l = min(jV /+l , Ny 12)). 

Ideally, the overall solution process (including the outer, time step iteration) should be defined using 
this reduced grid. However, as an intermediate step in testing the MG solution method, the approach 
taken has been to find some “reduced-grid” solution approximation for the full problem For this, the 
original grids are used, but it is desired to keep an approximation that lies in the corresponding reduced 
space. That is, the solution approximation is defined on all of G 1 , but has interpolated values on those 
points not in G 1 . The goal, perhaps somewhat artificial, has been to modify the MG algorithm so that 
the solution naturally remains in the reduced space while maintaining good MG convergence, as 
measured by the fine grid resid ual . That is, if the solution to the actual problem lies in this space, the 
MG method should converge to that solution with usual MG efficiency. To this end, the first step was 
to start with the full problem, and during plane relaxation near the poles, simply skip relaxation on 
levels finer than the new, reduced resolution there. That is, if >c/ > 1, then within the PMG algorithm, 

do not perform relaxation sweeps on levels m = 1, ... , log 2 // y . If the initial approximation on the plane 
is contained in this reduced solution space, then the corrected solution after the plane MG cycle will 
also be in this space. 

Using the above algorithm, one further requirement on the grid reduction is given by the global 
coarsening for practical reasons. Horizontal bilinear interpolation from level /+1 to level / should give a 
solution approximation that stays in the reduced space. That is, along each A-line in a particular /-plane, 
it should be a linear interpolant from N\ points along that line. If / is even, this is automatic, since the 

correction obtained from level /+1 is the interpolant from A^ /2 < N j points. If/ is odd, however, since 
plane / interpolates from planes // 2 and jf. 2 +1 (where //2 actually indicates the greatest integer iny/2), 
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the correction appears as a linear interpolant from max( N 1 *^ , A^'+i ) points. This requires that 
Nj > ma\(Np 2 , A^'+i ) • Taking the first value, this means that 

N 1 1 (Jj = N‘j > N';± = min(yV ,+1 ,/V'„,/2) 

= min( N'/2 , N l j 2/^_, ) 

Thus, f/j < 2 Similarly, < 2f/ J+] . This basically precludes large jumps in grid reduction, at least 
from an even to an odd plane. (The poles represent a “naturally” reduced plane with fJ N = f! s = N 1 , so 
there is no restriction in jumps from the poles to adjacent planes.) In the test cases here, a factor of 2 
reduction from one plane to the next was natural, so no problems were encountered. It should be noted 
that a factor of 2 from one plane to the next on level / corresponds to a factor of 4 on level /+1 , so for 
fine enough horizontal meshes, problems could be encountered and further restrictions would be 
necessary. However, in the meshes used and the reduction required here, the conditions above were 
met on all levels. 

It was noted that, in this approach, one goal was to maintain good MG convergence as measured by 
reduction of the fine grid residual. This is quite restrictive, and as noted, somewhat artificial. A more 
realistic method would be to simply measure the residual on the reduced grid problem, where the right- 
hand side is fixed. In the approach taken here, simply skipping relaxation on finer levels in plane 
relaxation results in recomputation of the right-hand sides of the “target” (coarser grid) plane problem. 
In testing convergence of the linear plane problems, it was found that, even when the finest grid right- 
hand side was set to zero (so that the solution was zero) and the initial error was constructed to lie in 
the reduced space, convergence per cycle degraded with the number of levels skipped in relaxation. 
Using a variational definition of the coarser grid linear operators as follows restored good convergence 
factors, as seen in Section 5.2. 

Z, m+1 = R” +l L m Pp for m = l, ..., M. (4.17) 

This has been implemented in the linear plane solvers only, and not in the global coarse grid 
discretizations. It is an open question whether the variational definition of the operators will actually be 
necessary when more realistic convergence measures for the reduced problem are introduced. 
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4.5 Efficiency Considerations. 


As the focus on work to date has been to obtain a working algorithm, concentrating on robustness 
rather than efficiency, timing comparisons for the method have not been included. Currently, the 
algorithm is probably much slower than comparable explicit or semi-implicit u-v models. In the 
algorithm presented here, the main computational work is in relaxation. Computation of the residuals is 
the second most expensive process. Some improvement can be gained here by more careful 
programming. However, the largest gains will be made in relaxation. This can be done in two ways. In 
both the line and plane relaxation, the main work is involved in the vertical line solves. The line problem 
is solved as a folly coupled system, so that the size of the matrix involved is approximately 6 K rows 
and columns. A direct solver is used that takes little advantage of the basic banded structure (although 
there is coupling at all layers to the surface pressure through the definition of <£). A more efficient solver 
can be written for this. The second method for making relaxation more efficient is to determine the real 
coupling of equations necessary. Often in MG applications, relaxation of the different unknowns can be 
done in separate sweeps. As a trivial example, consider the case of two grid functions (unknowns) 
defined on the mesh, in which each satisfies a Poisson equation. Since there are two independent 
problems defined, relaxation for each can be separate, but the type of relaxation grouping (line, plane, 
etc.) is determined by the geometry. More often, the connection between unknowns is of lower- 
triangular form, so that the unknowns can be relaxed in sequence, smoothing the error in each 
separately without adversely affecting residuals of previously relaxed equations. As noted, the current 
algorithm was written for robustness rather than efficiency, and more system analysis is required to 
determine the minimal amount of coupling necessary for relaxation (that is, which equations or groups 
of equations can be relaxed separately for the different quantities). Complete separation (which is 
unrealistic) would result in the direct solution of 6 KxK matrix problems per vertical line solve, rather 
than one 6Kx6K matrix problem, resulting in considerable savings. More realistically, it can be 
expected that some coupling will be necessary. However, any decoupling, combined with improved 
matrix solution routines tailored for the task, will improve efficiency greatly. 
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5. Numerical Results. 


In this section, numerical results are presented showing the behavior of the algorithm on several sample 
problems at single timesteps. Such problems are still being examined, and no results are presented 
showing the behavior of the algorithm in the M time-dependent model. In Section 5. 1 , the several 
model problems are presented, and the results with the MG algorithm are shown. There are problems 
for which the basic method of Section 4 does not converge. These are presented in Section 5.2. The 
difficulties with one of these problems is examined in detail, and a modified algorithm, using the 
reduced grids presented in Section 4.4, is shown to converge there. 

5.1 Initial Tests. 

The first set of problems considered are as follows: 

Problem 1. This example is a model system with no orography and a representative solution 
approximation. The mesh is defined by Nx = 64, N# = 32, and K= 18. A sigma vertical coordinate 
system (3.32) is used with uniform vertical mesh spacing (i.e., £ = k!K). 

Problem 2. This example is a model system with no orography and a representative solution 
approximation. The mesh is defined by N\ — 64, N# — 32, and K — 1 8. A hybrid vertical coordinate 
system (3.33) is used with a = 20 and 6U = 265. Different resolution is used in the lower and upper 
layers, with £ = 265, £ - £., = 6, for k = 1 , 2, ... 1 5, and £ - £_i = 40 otherwise, so that £»* = 475. 


Problem 3. This example is the same as the previous problem with the exception of the definition of the 
vertical coordinate. The hybrid coordinate (3.33) is still used, but with ct— 10 and dmn = 265. This 
time, £ = 265, £-£., = 8 2 / 3 for k= 1,2, ... 15, and £ - £.i = 26 2 / 3 otherwise, so that again £*x = 
475. 

First, the behavior of line and plane relaxation is examined (for Problems 1 and 2 only). The following 
tables show the behavior of the plane solvers used. Tables 1 -3 contain results for Problem 1 . For 
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j (plane) 

Eq. 1 

Eq. 2 

Eq. 3 

Eq. 4 

Eq. 5 

Eq. 6 

32 

0.380 

0.380 

0.382 

0.380 

0.384 

0.380 

31 

0.092 

0.089 

0.087 

0.092 

0.092 

0.085 

30 

0.089 

0.084 

0.078 

0.089 

0.085 

0.080 

25 

0.052 

0.052 

0.037 

0.052 

0.052 

0.038 

17 

0.030 

0.031 

0.018 

0.031 

0.031 

x r* tx 

0.017 

i- 1 1 


Table 1 . Average convergence factors for MG plane solver (linear problem) for Problem 1. 


j (plane) 

Eq. 1 

Eq. 2 

Eq.3 

Eq. 4 i 

Eq. 5 

Eq. 6 

32 

0.957 

0.979 

0.979 

0.957 

1.001 

0.966 

31 

0.966 

0.928 

0.928 

0.966 

0.947 

0.890 1 

30 

0.869 

0.850 

0.850 

0.869 

0.861 

0.774 

25 

0.445 

0.445 

0.445 

0.445 

0.447 

0.299 

17 

0.241 

0.250 

0.250 

0.250 

0.250 

0.146 


Table 2. Average convergence factors for relaxation only (linear problem) for Problem 1 . 


j (plane) 

Eq. 1 

Eq. 2 

Eq.3 

Eq. 4 

Eq. 5 

Eq. 6 

32 

0.377 

0.263 

0.891 

0.377 

1.000 

1.006 

31 

0.091 

0.090* 

0.994 

0.084 

1.012 

1.001 

30 

0.080 

0.083* 

1.003 

0.083 

1.008 

0.994 

25 

0.053* 

0.052* 

1.005 

0.050* 

0.996 

0.989 

17 

0.031* 

0.031* 

1.000 

0.031* 

0.992 

0.975 


Table 3. Average convergence factors for MG plane solver (nonlinear problem) for Problem 1 . 


Table 


j (plane) 

Eq. 1 

Eq. 2 

Eq.3 

Eq. 4 

Eq. 5 

Eq. 6 

32 

0.083 

0.061 

0.079 

0.081 

0.082 

0.083 

31 

0.082 

0.089 

0.083 

0.072 

0.087 1 

0.090 

30 

0.074 

0.084 

0.082 

0.086 

0.079 

0.081 

25 

0.055 

0.052 

0.035 

0.055 

0.051 

0.052 

17 

0.031 

0.031 

0.017 

0.032 

0.032 

• i \ n- tv 

0.031 

it 


4 . Average convergence factors for MG plane solver (linear problem) for Problem 2. 


j (plane) 

Eq. 1 

Eq. 2 

Eq.3 

Eq. 4 

Eq. 5 

Eq. 6 

32 

0.956 

0.979 

0.965 

0.925 

0.972 

0.967 

31 

0.978 

0.928 

0.877 

0.970 

0.948 

0.955 

30 

0.870 

0.850 

0.755 

0.865 

0.865 

0.862 

25 

0.445 

0.445 ^ 

0.285 

0.445 

0.445 

0.440 

17 

0.250 

0.250 

0.141 

0.250 

0.251 

0.251 


Table 5. Average convergence factors for relaxation only (linear problem) for Problem 2 
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j (plane) 

Eq. 1 

Eq. 2 

■301 

Eq. 4 


Eq. 6 

32 

0.085 

0.191 

1.053 

0.084 

0.979 

0.812 

31 

0.087 

0.130 

1.009 

0.083 

1.049 

0.000 

30 

0.078 

0.084 

1.025 

0.081 



25 

0.052 

0.055* 

0.739 

0.036* 

1.075 

0.909 

17 

0.032* 

0.033* 

0.984 

0.017* 

1.055 

0.922 


Table 6. Average convergence factors for MG plane solver (nonlinear problem) for Problem 2. 


V-cycle results for the full problems are given in Tables 7, 8, and 9, for problems 1, 2, and 3, 
respectively. These contain initial residuals for each equation and residuals after each cycle, and an 
effective convergence factor (the maximum convergence factor for equations 1, 2, and 4). For reasons 
discussed below, the norm of the right-hand side of each equation is also given. In both cases, 
convergence factors appear satisfactory for a few cycles, then the residuals stagnate. This type of 
behavior can often be attributed to one (or a combination) of three causes: The level of roundoff error 
can be reached; the solver may not converge for certain error components; or the right-hand side may 
not be in the range of the operator. The level of roundoff error has clearly not been reached, since that 
should ideally give residuals of a size 10' 16 -|j/||. Here, the third appears the most likely. It is clear that 
the operator itself is singular, since there are no horizontal boundary conditions and all % and ¥ only 
occur in conjunction with horizontal derivatives. Thus, arbitrary constants can be added to ^and y/ at 
each vertical level without affecting the residual. This means that the range of the operator cannot be 
the set of all possible right-hand sides, and some consistency condition must be met if the problem is to 
have a solution. For these test problems, all grid functions used (the right-hand sides, the initial solution 
approximation, and various coefficients) were given to six decimal places only, so it seems reasonable 
to assume that the residual could be reduced to approximately 10 _6 *||/|| without consistency problems. 
Comparison of the largest residual (equation 4) with ||/ 4 1| shows that this appears to be the case. (For 
smaller residuals, particularly equation 2, it is likely that lack of convergence of the other equations 
could interfere there, so that reduction of those residuals do not reach the expected levels.) 


Unfortunately, for a complicated nonlinear system such as this, it is difficult to determine the proper 
consistency condition and to somehow enforce it. As a test, 20 cycles were performed, and a new 
right-hand side computed consistent with the solution obtained (that is, the right-hand side was simply 
set to the operator applied to the solution). Then, the original initial approximation was restored and 
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Cycle 

Eq. 1 

Eq. 2 

Eq. 3 

Eq. 4 

Eq. 5 

Eq. 6 

Factor 

ll/ll 

O.OOE+OO 

1.52E-09 

O.OOE+OO 

4.76E+-02 

4.76E-02 

O.OOE+OO 

- 

Initial 

4.03E-10 

8.75E-10 

4.16E-02 

3.60E-01 

3.04E-04 

6.88E-09 

- 

1 

9.19E-1 1 

1.16E-10 

1.71E-05 ] 

3.17E-02 

3.18E-18 

1.96E-16 

0.088 

2 

1.14E-1 1 

1.08E-1 1 

1.33E-08 

1.90E-03 

3.22E-18 

2.43E-19 

0.060 

3 

1.71E-12 * 

1.51E-12 

1.22E-1 1 

2.73E-04 

3.19E-18 

2.41E-19 

0.144 

4 

2.76E-13 

7.71E-13 

2.22E-12 l 

1.27E-04 

3.20E-18 

2.39E-19 

0.465 

5 

6.99E-14 

7.57E-1 3 1 

2.23E-12 

1.23E-0<P 

3.20E-18 

2.41E-19 

0.969 

6 

4.93E-14 

7.57E-13 

2.22E-12 

1.23E-04 

3.19E-18 

2.44E-19 

1.000 

7 

4.83E-14 

7.57E-13 

2.23E-12 

1.23E-04 

3.21E-18 

2.40E-19 

1.000 

8 

4.83E-14 

7.57E-13 

2.22E-12 

1.23E-04 

3.17E-18 

2.44E-19 

1.000 

9 

4.83E-14 

7.57E-13 

2.22E-12 

1.23E-04 

3.21E-18 

2.42E-19 

1.000 

10 

4.83 E- 14 

7.57E-13 

2.21E-12 

1.23E-04 

3.23E-18 

2.41E-19 

1.000 


Table 7. Convergence history for Problem 1 (original forcing terms). 


Cycle 

Eq. 1 

Eq. 2 

Eq. 3 

Eq. 4 

Eq -5 

Eq. 6 

Factor 

ii/ii 

0.00E+00 

9.12E-10 

0.00E+00 

5.74E+00 

3.30E-02 

O.OOE+OO 

- 

Initial 

1.21E-1 1 

1.96E-09 

3.27E-04 

5.65E-05 

6.73E-09 

2.09E-03 


i 

1.71E-09 

2.39E-10 

1.89E-06 

1.96E-03 

2.15E-18 

2.82E-08 

>1 

2 

5.94E-10 

7.58E-1 1 

9.09E-09 

1.95E-04 

1.99E-18 

1.40E-12 

0.347 

3 

1.91E-10 

2.64E-1 1 

3.24E-1 1 

4.71E-05 

1.89E-18 

3.01 E- 1 4 

0.348 

4 

6.07E-1 1 

1.1 IE-1 1 

3.87E-13 

2.46E-05 

1.94E-18 

2.51E-15 

0.522 

5 

1.92E-11 

8.32E-12 

1.17E-13 

1.95E-05 

1.93E-18 

8.76E-16 

0.793 

6 

6.03E-lz 

8.10E-12 

6.47E-14 

1.86E-05 

1.89E-18 

5.86E-16 

0.954 

7 

1.90E-12 

8.11E-12 

5.26E-14 

1.83E-05 

1.94E-18 

5.23E-16 

0.984 

8 

6.25E-13 

8.12E-12 

4.92E-14 

1.83E-05 

1.91E-18 

5.00E-16 

1.000 

9 

2.70E-13 

8.12E-12 

4.78E-14 

1.83E-05 

1.95E-18 

4.99E-16 

1.000 

10 

1.98E-13 

8.12E-12 

4.79E-14 

1.82E-05 

1.95E-18 

4.97E-16 

0.995 


Table 8. Convergence history for Problem 2 (original forcing terms). 


Cycle 

Eq. 1 

Eq. 2 

Eq. 3 

Eq. 4 

Eq. 5 

Eq- 6 

Factor 

II f\\ 

O.OOE+OO 

6.48E-10 

O.OOE+OO 

4.62E+00 

2.85E-02 

O.OOE+OO 

- 

Initial 

2.30E-10 

1.14E-09 

2.91E-04 

1.90E-03 

1.38E-05 

3.57E-04 

- 

i 

4.70E-1 1 

1.01E-10 

1.11E-07 

3.43E-04 

1.54E-18 

1.12E-10 

0.204 

2 

9.55E-12 

8.39E-12 

2.78E-10 

1.96E-05 

1.51E-18 

6.02E-14 

0.203 

3 

2.73E-12 

1.37E-12 

9.57E-13 

2.85E-06 

1.54E-18 

3.86E-16 

0.286 

4 

7.62E-13 

1.06E-12 1 

1.37E-14 

1.61E-06 

1.52E-18 

1.58E-16 

0.774 

5 

2.06E-13 

1.04E-12 

1.32E-14 

1.68E-06 

1.52E-18 

1.58E-16 

1.043 

6 

5.53E-14 

1.04E-12 

1.33E-14 

1.69E-06 

1.53E-18 1 

1.60E-16 

1.006 

7 

2.34E-14 

1.04E-12 

1.32E-14 

1.69E-06 

1.54E-18 

1.56E-16 

1.000 

8 

2.20E-14 

1.04E-12 

1.34E-14 

1.69E-06 

1.54E-18 

1.59E-16 

1.000 

9 

2.27E-14 

1.04E-12 

1.33E-14 

1.69E-06 

1.53E-18 

1.60E-16 

1.032 

10 

2.29E-14 

1.04E-12 

1.32E-14 

1.69E-06 

1.51E-18 

1.56E-16 

1.009 


Table 9. Convergence history for Problem 3 (original forcing terms). 
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the problem re-solved. Results for problems 1, 2, and 3 are shown in Tables 10, 11, and 12, 
respectively. Convergence now appears quite good. This cannot be viewed as conclusive, however, 
since this test may be construed as cheating in a way. If there are error components that are untouched 
by the MG solution process, it could be simply that the new problem is manufactured so that these 
error components are not present. However, this seems unlikely. This type of behavior is consistent 
with that for simpler linear systems in which the operator is singular, such as the Poisson problem with 
foil Neumann boundary conditions. There, the same type of singularity exists, and in an inconsistent 
problem, MG will converge well to a solution to the consistent problem, although as measured by the 
residual, convergence stagnates. The fact that the modified problem continues to converge well long 
after the original solution process appears to stagnate indicates that something similar is occurring here, 
and stagnation is due to inconsistency of the right-hand side. 


Problem 1 and 3 convergence histories generally look reasonable, while problem 2 residuals show a 
jump at the first cycle. This may indicate that the solution process actually converges to a different 
solution than the desired one. The main difference between problems 2 and 3 is the definition of a used 
in the hybrid coordinate. This needs closer examination, and awaits further study. 


Cycle 

Eq. 1 

Eq. 2 

Eq. 3 

Eq. 4 

Eq. 5 

Eq. 6 

Factor 

H/ll 

0.00E+00 

1.52E-09 

O.OOE+OO 

4.76E+02 

4.76E-02 

0.00E+00 

- 

Initial 

4.03E-10 

8.75E-10 

4.16E-02 

3.60E-01 

3.04E-04 

6.88E-09 

- 

1 

9.19E-1 1 

1.16E-10 

1.71E-05 

3.16E-02 

3.21E-18 

3.2 IE- 18 | 

0.124 

2 

1.14E-1 1 

1.08E-11 

1.33E-08 

1.90E-03 

3.19E-18 

3.19E-18 ! 

0.150 ^ 

3 

1.71E-12 

1.33E-12 

1.22E-1 1 

2.43E-04 

3.22E-18 

— 

3.22E-18 

0.159 

4 

2.72E-13 

1.66E-13 

2.24E-12 

3.00E-05 

3.21E-18 

3.21E-18 

0.192 

5 

5.21E-14 

2.22E-14 

2.24E-12 

4.01E-06 

3.21E-18 

3.21E-18 

0.228 

6 

1.19E-14 

3.19E-15 

2.23E-12 

5.59E-07 

3.22E-18 

3.22E-18 

0.250 

7 

2.98E-15 

5.53E-16 

2.22E-12 

8.82E-08 

3.23E-18 1 

3.23E-18 

0.258 

8 

7.70E-16 

1.19E-16~1 

2.23E-12 

1.70E-08 

3.23E-18 

3.23E-18 1 

0.261 

9 

2.01E-16 

2.92E-17 

2.22E-12 

I 3.89E-09 

3.21E-1© 

3.21E-18 

0.262 

10 

5.26E-17 

7.50E-18 

2.21E-12 

9.75E-10 

3.20E-18 

3.20E-18 

0.262 

11 

1.38E-17 

1.95E-18 

2.24E-12 

2.52E-10 

3.18E-18 

3.18E-18 

0.262 

12 

3.61E-18 

5.12E-19 

2.11E-12 

6.58E-1 1 

3.20E-18 

3.20E-18 

0.262 

13 

9.46E-19 

1.37E-19 

2.00E-12 

1.72E-11 

3.19E-18 

3.19E-18 

0.265 

14 

2.48E-19 

4.42E-20 

1.82E-12 

4.55E-12 

3.23E-18 

3.23E-18 

0.629 

15 

6.49E-20 

2.78E-20 

1.57E-12 

1.28E-12 

3.19E-18 

3.19E-18 

0.942 


Table 10. Convergence history for Problem 1 (modified forcing terms). 
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0.00E+00 

9.12E-10 

O.OOE+OO 

1 .21 E-l 1 

1.96E-09 

3.27E-04 

1.71E-09 

2.39E-10 

1.89E-06 

5.94E-10 

7.60E-1 1 

9.08E-09 

1.91E-10 

2.58E-1 1 

3.24E-1 1 

6.07E-1 1 

8.31 E-l 2 

1.98E-13 

1.92E-1 1 

2.63E-12 

2.79E-14 

6.01E-12 

8.24E-13 

1.95E-14 

1.87E-12 

2.55E-13 

1.92E-14 

5.76E-13 

7.85E-14 

1.90E-14 


2.39E-14 

1.91E-14 

5.34E-14 

7.22E-15 

1.90E-14 

1.60E-14 

2.16E-15 

1.88E-14 

4.78E-15 

6.42E-16 

1.90E-14 

1.41E-15 

1.89E-16 

1.89E-14 

4.14E-16 

5.52E-17 

1.90E-14 

1.20E-16 

1 1.60E-17 

1.93 E-l 4 


E 


5.74E+00 


5.93E-05 


1.97E-03 


1.87E-04 


4.06E-05 


1.30E-05 


4.10E-06 


1.29E-06 


4.03 E-07 


1.25E-07 


3.84E-08 


1.17E-08 


3.53E-09 


3.13E-10 


9.19E-11 


2.68E-1 1 


E 


3.30E-02 


6.73E-09 


2.16E-18 


1.94E-18 


1.93E-18 


1.91E-18 


1.96E-18 


1.89E-18 


1.95E-18 


1.95E-18 


1.97E-18 


2.00E-18 


1.97E-18 


1.95E-18 


E 


0.00E+00 


2.09E-03 


2.82E-08 


1.40E-12 


1.23E-15 


2.02E-16 


1.50E-16 


1.48E-16 


1.43E-16 


1.39E-16 


1.24E-16 


1.17E-16 


1.14E-16 


1.09E-16 


1.04E-16 


Factor 


>1 


0.095 


0.217 


0.320 


0.315 


0.315 


0.312 


0.310 


0.307 


0.305 


0.302 


0.297 


0.298 


0.294 


0.292 


Table 11. Convergence history for Problem 2 (modified forcing terms). 



2.30E-10 


4.70E-1 1 


9.56E-12 


2.74E-12 


7.70E-13 


2.13E-13 


5.83E-14 


1.59E-14 


4.35E-15 


1.19E-15 


3.24E-16 


6.56E-18 


1.79E-18 


4.86E-19 


E 


6.48E-10 


1.14E-09 


1.01E-10 


8.35E-12 


8.86E-13 


1.29E-13 


2.91E-14 


7.87E-15 


2.17E-15 


5.95E-16 


1.63E-16 


4.45E-17 


1.21E-17 


3.31E-18 


9.01 E-l 9 


2.47E-19 


7.20E-20 


E 


1.33E-14 


2.91E-04 


1.11 E-07 


2.78E-10 


9.57E-13 


1.37E-14 


1.33E-14 


1.33E-14 


1.32E-14 


1.32E-14 


1.32E-14 


1.33E-14 


1.32E-14 


1.32E-14 


1.28E-14 


1.25E-14 


1.23E-14 


E 


4.62E+00 


1.90E-03 


3.41 E-04 


2.07E-05 


1.88E-06 


2.17E-07 


4.96E-08 


1.27E-08 


3.45E-09 


9.38E-10 


2.56E-10 


1.90E-1 1 


5.17E-12 


1.41E-12 


3.84E-13 


1.04E-13 


Eq.5 


2.85E-02 


1.38E-05 


1.56E-18 


1.54E-18 


1.53E-18 


1.51E-18 


1.52E-18 


1.53E-18 


1.54E-18 


1.53E-18 


1.50E-18 


1.53E-18 


1.53 E-l 8 


1.52E-18 


1.55E-18 


1.56E-18 


E 


0.00E+00 


3.57E-04 


1.12E-10 


6.02E-14 


4.12E-16 


1.59E-16 


1.58E-16 


1.60E-16 


1.56E-16 


1.49E-16 


1.40E-16 


1.23E-16 


1.09E-16 


9.72E-17 


8.66E-17 


7.67E-17 


7.03E-17 


Factor 



0.204 


0.203 


0.287 


0.281 


0.277 


0.274 


0.276 


0.274 


0.274 


0.273 


0.273 


0.274 


0.273 


0.274 


0.291 


Table 12. Convergence history for Problem 3 (modified forcing terms). 


39 



































































































































































































































































5.2 Further Tests and the Reduced grid Problem 


During the development of the MG solver, difficulties with convergence led to generating test 
problems in which the atmosphere was in a state of rest (that is, u = v = £= 0). Here, the solution 
should remain unchanged, and the purpose of these problems was to make sure that the MG rolver 
maintained this state of rest. Problems arose with the more highly refined mesh, . The problem, and the 
approach taken to deal with it, are discussed in this section. 

Problem 4. This example is a model system with no orography. The mesh is defined by M, ■ = 64, N t = 

32, and K - 1 8. A sigma vertical coordinate system ( 3 . 32 ) is used with uniform vertical mesh spacing 

(i.e., & = kl K). 

Problem 5. This example is a model system with no orography. The mesh is defined by N> = 64, N, - 
32, and K = 18. Ahybrid vertical coordinate system (3.33) is used with 5 and A. = 260. Uniform 

vertical mesh spacing is used with <§> = 260 and = 550. 

Problem 6. This example is a model system with no orography. Here, a more highly refined horizontal 
mesh is used, with with% - 128, JV,= 64, and K= 18. A sigma vertical coordinate system (3.32) is 

used with uniform vertical mesh spacing (i.e., £ - klK). 

For problems 4 and 5, the solver did not diverge, maintaining the state of rest. However, The solver 
Med on Problem 6, when a negative value of pressure was obtained (which makes the computatton of 
n impossible). It was found that the plane MG solver was diverging near the poles. This is 
demonstrated Table 1 3, which shows the behavior of plane relaxation on the linear problem for) - 63 
(adjacent to the north pole), 62, 61, 60, 59, 48 (45' N), and 32 (the equator). 
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j (plane) 

HI 

esh 

ECU 

KSBfl 

esh 

ESEfli 

63 

120.3 

39.43 

122.2 

120.3 

123.2 

121.1 

62 

21.50 

0.945 

20.74 

21.51 

21.52 

4.989 

61 

1.667 

0.947 

1.642 

1.667 

1.555 

0.924 

60 

0.839 

0.923 

0.881 

0.838 

0.977 

0.904 

59 

0.891 

0.889 

0.853 

0.891 

0.935 

0.859 

48 

0.478 

0.444 

0.404 

0.478 

0.456 

0.397 

32 

0.250 

0.250 

0.221 

0.269 

0.258 

0.218 


Table 13. Average convergence factors for line relaxation only, Problem 6, full linear problem. 


Trial and error was used to determine the equations and terms responsible for this (that is, to identify 
the simplest subsystem that still exhibiting this divergent behavior in the linear solver). The resulting 
system retained Equations (3.13), (3.14b), (3.15), and (3.16), and the unknowns {0\ M). Note 

that M 0 , and (through the boundary condition (3.30)) p 0 , have been excluded. The resulting equations, 
showing the terms retained, are as follows. 


Equation (3.13) for h=\, 2, ..., 

K: 





+ v 2 ff ] r / t 2 

(5.1) 

Equation (3.14b) for k= 1, 2, .. 

,K- 1 : 






(5.2) 

Equation (3.15) for k= 1, 2, ..., 

K: 





(v 7 z), =/,* 

(5-3) 

Equation (3.15) for k=0, 1 , 2, 

...JC: 






(5.4) 


Here, n is frozen using the current values of p, all values of rare equal, and vertical mesh spacing is 
uniform (so that solid and virtual levels coincide) . Now, repeating the plane relaxation tests on this 
reduced system gives the results shown in Table 14. Note that both diverge for / = 61, 62, and 63, 
while they converge elsewhere. 
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j (plane) 

Eq. 1 

Eq. 2 

Eq. 3 

Eq. 4 

Eq. 5 

Eq. 6 

63 

- 

21.09 

21.03 

21.07 

21.06 

- 

62 1 

- 

44.73 

44.71 

44.77 

44.77 

- 

61 

- 

1.440 

2.335 

2.339 

2.335 

- 

60 

- 

0.928 

0.929 

0.929 

0.926 

- 

59 

i 

0.892 

0.895 

0.891 

0.893 

- 

48 


0.445 

0.441 

0.443 

0.441 

- 

32 

- 

0.250 

0.250 

0.249 1 

0.250 

- 


Table 14. Average convergence factors for line relaxation only. Problem 6, reduced linear problem. 


To see why divergence occurs, this system is first simplified as follows: 

1 . Eliminate % in Equation (5.1) using (5.3). 

2. Eliminate 6 in Equations (5.1) and (5.2) using Equation (5.4). 

3. Apply 8c to Equation (5.1), r V 2 to Equation (5.2), and subtract to eliminate M. 
Explicitly writing the vertical differences and averages, the resulting equation in £ has the form: 


§k±i ~2<£t + §k-i . j 

hi 2 


n*+i n* a y2 e . n* n*-i p y 2 ? 

] l k+ 1 V SA+1+ 7 A *-1 V bk - 1 

he hf. 


= f (5.5) 


‘i '*£ 

for A: = 1, 2, ..., AM. The boundary conditions are as before, with £p = %k = 0. Now, restricted to a bi- 
plane j, the operator V 2 has the form 

+c ' A -^ +c ‘ A >-'^ ’ (5 - 6) 


where 


c w = c e = ■ 

C 7 


r 2 A 2 


cos 2 <p. 


and 


(r w 1 c e ) — * J+V2 + COS ^ lU2 
1 1 J a 2 h 2 vos <f>j 


Thus, the operator can be scaled and written in stencil form as: 


'0 

-1 

O' 


'-i 

2 + s -\ 


0 

0 

O' 

0 

2 

0 

+4 

0 

0 

0 

+ C k 

0 

0 

0 

0 

-1 

0_ 


0 

0 

0 


-1 

2 + e -lj 


where 

* 2 a 2 h 2 cos 2 hg 

f t _, n,-ry, 

* 2a 2 /i 2 cos 2 <j>j ht 


(5.7) 


(5.8) 


(5.9) 


and 
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£=cos</>j(cos</> j+V2 +COS ^_ 1/2 ). (5.10) 

As before, h = n IN. Also, note that at grid lines adjacent to the poles (/' = 1 and AM), cos$ » h, so that 
2 h 2 . Tests show that line relaxation on a model system (5.5) with c* = c\ = c for all k and e= 0 
diverges for c > 0.5. In particular, note that, for this simplified model, if c = 0.5, the resulting stencil is: 

'-'2 0 -V 
020 
_-* 2 0 - f 2 

This is the so-called skewed Laplace operator, which has an MG convergence fector of 1 .0 when 
semicoarsening (or standard coarsening, lor that matter) is used, since there is no smoothing in the left- 
right direction. In practice, there is some convergence due to relaxation because of the Dirichlet 
boundary conditions at k = 0 and k = K. For acceptable convergence, the value should actually be 
somewhat less than 0.5. In Problem 6 , the c*’s increase with £ and decrease according to the square of 
the distance from the pole. The maximum value of c k is 7.05 at 7 = 63, 1 .76 at 7 = 62, 0.78 at/ = 61, 

0.31 at 7 = 60, and 0.29 at 7 =59. Thus, divergence for this simplified model corresponds to that for the 
plane relaxation on the reduced equations and the full system. It should be noted that results are 
symmetrical about the equator, so the same behavior is observed at/ = 1 as at 7 = 63, for example. 

It appears that the stencil is some discretization of 3 •/ +cV 2 , which is nicely elliptic, but that the 
vertical positioning of the equations and unknowns in the discretization leads to a vertical averaging for 
the cV 2 term that introduces non-ellipticity near the poles. Several methods were tested for dealing 
with this problem, all based on the fact that the grid is basically over-resolved in the ^-direction near 
the poles. Various A-direction averaging methods were tested, both averaging the solution to reduce 
high frequency error and introducing using averaging in constructing the horizontal differences. These 
were unsuccessful, however. Another method was then used in which longer-range larger stencils were 
used for the difference operators in the /.-direction, which is described in Section 4.4. 

Note that the denominators in the first terms in (5.9) correspond to twice the square of the horizontal 
mesh size on the spherical mesh, so that when the reduced grid method is used, the ct s for plane 7 are 
are reduced by a fector of ff r Using the values for c* obtained for Problem 6 , minimum values of/} 
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needed for convergence are using p - p,i - 4, pi - psi - 2, and p 1 otherwise. Convergence of 
plane MG relaxation for j = 63, 59 for values of/} = 1, 2, 4, and 8 are given in Table 15 for the 
reduced system (5.1)-(5.4), and for the full system (3.12)-(3.17) in Table 16. Here, standard coarse 
grid discretizations for the operators are used. Note that divergence does occur where predicted. 
However, as noted in Section 4.4, there appears to be a limit to convergence factors that degrades as p 
is increased. The convergence factors are measured on the finest level, not the actual reduced grid 
level, so this could be an artifact of this. Nevertheless, it was found that the use of variational operators 
(4.17) improved convergence, avoiding this limitation. Corresponding results for the variational 
operators are shown in Table 1 7 for the reduced system (5.1)-(5.4), and in Table 1 8 for the full system 
(3.12X3.17). 


j (plane) 

Pi~ 1 

p = 2 

A = 4 

A = 8 

63 

1474. 

1502. 

0.600 

0.630 

62 


0.483 

0.600 

0.630 

61 

4.445 

0.483 * 

0.600 

0.630 

■K3H ■ 

0.104 

0.483 

0.600 

0.630 

59 

0.095 

0.483 

0.600 

0.630 


Table 15. Overall average convergence factors for plane MG cycles for Problem 6, 
(5.1)-(5.4), standard CG operators, using the reduced grids defined by p. 


reduced system 


j (plane) 

A = 1 

A = 2 

A = 4 

A~ 8 

63 

8096. 

301.0 

0.599 

0.626 

62 

222.7 

0.481 

0.601 

0.626 

61 1 

2.508 

0.480 

0.600 

0.626 

60 

0.101 

0.481 

0.600 

0.626 

59 

0.096 

0.481 

0.600 

0.626 


Table 16. Overall average convergence factors for plane MG cycles for Problem 6, 
(3.17), standard CG operators, using the reduced grids defined by p. 


full system (3.12)- 


j (plane) 

A= 1 

A 2 

A = 4 

A = 8 _ 

63 

>9999 

6253. 

8.796 

0.087 

62 


3.309 

0.077 

0.087 

61 

4.452 

0.124 

0.069 

0.087 


0.106 

0.077 

0.069 

0.087 

59 

0.087 

0.068 

0.069 

0.087 


(5.1)-(5.4), variational CG operators, using the reduced grids defined by p. 
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j (plane) 

P>~ 1 

P, = 2 

£lZA 

Pi § 

63 

>9999 

>9999 

6.168 

0.087 

62 

271.3 

2.437 

0.077 

0.087 

61 

2.490 

0.104 

0.069 

0.087 1 


0.094 

0.077 

0.069 

0.089 

59 

0.090 

0.068 

0.069 

0.087 


Table 18. Overall average convergence factors for plane MG cycles for Problem 6, 
(3.17), variational CG operators, using the reduced grids defined by p. 


full system (3.12)- 


One thing to note is that, with the variational discretization, the plane problems diverge for pti 4 and 
— 2 5 while convergence was obtained with standard discretizations. That is due to the feet that, for 
the vertical derivatives, the variational definition uses a horizontal averaging over three vertical grid 
lines, which affects smoothing and places somewhat stricter requirements the c*’s. Thus, A 3 = 8 and 
A 2 = 4, are required for convergence of the plane MG solver. 


Thus, the reduced grid method was implemented using for the fell Problem 6 using values p\-pbi~%, 
p^ = p ^2 = 4, a ~ A \ = 2, and p, = 1 otherwise. Convergence factors for the first five cycles are shown 
in the Table 19 below. The divergence problem has been solved, and it appears that the solution as 
given is accurate to close to the level of roundoff error. 


Cycle 

Eq. 1 

Eq. 2 

Eq. 3 

Eq. 4 

Eq. 5 

Eq. 6 

II /'ll 

0.00E+00 

0.00E+00 

0.00E+00 

2.49E+02 

O.OOE+OO 

O.OOE+OO 

Initial 

5. 69 E- 17 

4.66E-18 

4.57E-08 

1.69E-09 

5.58E-27 

7.91E-19 1 

1 

1.57E-17 

1.60E-18 

1.37E-10 

5.47E-09 

3.24E-27 

1.18E-19 

2 

7.15E-18 

7.33E-19 

9.34E-1 1 

3.45E-09 

3.44E-28 

1.26E-19 

3 

5.91E-18 

7.42E-19~1 

8.81E-1 1 

3.39E-09 

3.65E-28 

1.34E-19 

4 

5.48E-18 

7.92E-19 

8.31E-1 1 

3.36E-09 

3.84E-28 

1.41E-19 

5 

5.28E-18 

8.42E-19 

7.84E-1 1 

3.34E-09 

4.11E-28 

1.47E-19 


Table 19. Convergence history for Problem 6. 


A note of warning about the general use of reduced grids is needed here. For fixed r and h f , it is 
possible to reduce h enough so that the desired p, is larger than N. This may not be a practical problem, 
since higher accuracy will usually require both increased vertical and horizontal resolution, as well as 

decreased timesteps. 
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Finally, it should be noted that the method Med on a test problem with orography, a sigma vertical 
coordinate system (3.32) with uniform vertical mesh spacing (i.e., & = kIK), and the mesh is defined by 
= 64, JV # = 32, and K = 1 8. The reason for this is still under study. 

6. Conclusions. 

An FAS FMG solver has been developed which shows promise for the atmospheric general circulation 
model based on the semi-Lagrangian advection of potential vorticity (PV) presented here. Additional 
development of the solver is needed, both for robustness and efficiency reasons, as there are test 
problems for which the method was unable to converge. An important remaining step is to test the 
solver in the full model in order. Since £-line relaxation dominates the computation, both in the global 
sweeps away from the poles and within the MG plane solvers, this is where the most speedup is to be 
gained. Further system analysis to determine the needed coupling of equations necessary in relaxation 
and the development of tailored linear solvers are expected to improve efficiency. 
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